function [err, VCresult]=VC_TC( sf, f)
% Analyze VC current time course (assume ONE pulse amplitude...)
%
global VOLTAGE CURRENT DFILE
global CONTROL

err = 0;

QueMessage('VC_TC: Preparing for analysis', 1);
dat = [];
time = [];
ivresult = []; % initialize it

% abstract general information
protocol=deblank(lower(CONTROL(sf).protocol));
rate = DFILE.rate*DFILE.nr_channel/1000;

[records,pts]=size(VOLTAGE);
% always compute the time of each stimulus, in seconds from the start of the data
a = DFILE.ztime;
TM = (a - a(1))/60; % time at which data was collected... gathered from ztime...

% Get analysis windows 
%stim_time=CONTROL(sf).stim_time; % get the stimulus array - usually 2 times....

% get the psp definition array. For field potentials, we have the following
% sets of times: (layer 2 DCN pf FP)
% n1b, n1e
% n2b, n2e
% p3b, p3e
% p4b, p4e
% It is assumed that the same response is repeated with similar latencies for each stimulus...


ts1=number_arg(CONTROL(sf).durho); % get the baseline duration (most of time before stimulus)

tsbs = number_arg(CONTROL(sf).durho)+number_arg(CONTROL(sf).deadwin);
tsb=number_arg(CONTROL(sf).durho)+number_arg(CONTROL(sf).durs1); % set window where most of analysis takes place)
tse=tsb+number_arg(CONTROL(sf).durs2);

% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
k = find(DFILE.rate < 1);
DFILE.rate(k) = 50;
RATES = (DFILE.rate * DFILE.nr_channel)/ 1000; % array of sampling rates, convert to msec

vc_base = floor((ts1 ./ RATES) - 0.5)
	vc_m1b=floor(tsbs./RATES); % for response
	vc_m1e=floor(0.5*tse./RATES);
   vc_m2b=floor(0.80*tse./RATES); 
   vc_m2e=floor(0.99*tse./RATES);

% measure holding current and "rmp"
for i = 1:records
	hold_cur(i) = mean(CURRENT(i,1:vc_base(i))); % save in array for later usage
	rmp_volt(i) = mean(VOLTAGE(i,1:vc_base(i)));
end;
CONTROL(sf).iHold = mean(hold_cur);
CONTROL(sf).Rmp = mean(rmp_volt);

% but first smooth the current out a bit
QueMessage('  --  Filtering');
for i = 1:records
	fsamp = 1000/RATES(i); % get sampling frequency
	fco = 200;		% cutoff frequency in Hz
	wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
	if(wco < 1) % if wco is > 1 then this is not a filter!
	   [b, a] = fir1(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
	   current(i,:) = filtfilt(b, a, CURRENT(i,:)')'; % filter all the traces...
   else
      current = CURRENT;
   end;
   
end
QueMessage('  --  Baseline subtraction');
for i = 1:records
	base = mean(current(i,1:vc_base(i))');
	current(i,:) = current(i,:) - base;
end;
QueMessage('  --  Measuring Currents');
% find min voltage in window 1 (N2, pop spike)
MPEAK_VC = zeros(records, 1);
MMEAN_VC = zeros(records, 1);

for i = 1:records
	MPEAK_VC(i) = max(current(i,vc_m1b(i):vc_m1e(i))'); % current peak 	flipped for direct comparison...
   MMEAN_VC(i) = mean(current(i,vc_m2b(i):vc_m2e(i))'); % current mean
end;

VC_TC.max=MPEAK_VC;
VC_TC.mean = MMEAN_VC;
VC_TC.time = TM;

% do F test on the windows...
%tbs=1;
%a=find(TM<600);
%tbe=a(length(a)); % use the last point accessible
%a=find(TM>900);
%tps=a(1);
%a=find(TM>1500);
%tpe=a(1);
%tbs
%tbe
%tps
%tpe
%size(TM)

%[P, F, dfR, dfF, yp1, yp2, xf, yf] = bears_f_test(TM(tbs:tbe)', MAX_FP2V(tbs:tbe), ...
%   TM(tps:tpe)', MAX_FP2V(tps:tpe));
%disp(sprintf('base: %7.2f mV  Sig: %7.2f mV, P=%6.3f', mean(MAX_FP2V(tbs:tbe)), ...
%   mean(MAX_FP2V(tps:tpe)), P))

CONTROL(sf).VC_TC = VC_TC;

QueMessage('  --  Analysis complete');


% plot if figure is set
if f>0
   h = findobj('Tag', 'VC_TC'); % check for pre-existing window
   if(isempty(h)) % if none, make one
      h = figure('Tag', 'VC_TC', 'Name', 'VC Time Course Analysis', 'NumberTitle', 'off');
   end
   figure(h); % otherwise, select it
   clf; % always clear the window...
   fsize = 7;
   msize = 3;
   tmax = max(time);    
   
   %plot Min voltage for N2
   subplot('Position',[0.1,0.30,0.8,0.650]);
   plot(TM, MPEAK_VC, '-ks', 'Markersize', msize, 'Markerfacecolor', 'k') % data in black
   hold on;
   plot(TM, MMEAN_VC, '-rx', 'Markersize', msize, 'Markerfacecolor', 'r'); % data in blue
   set(gca, 'FontSize', fsize);
   ylabel('VC current (nA)');
   grid;
   u=get(gca, 'YLim');
   set(gca, 'Ylim', [0 u(2)]);
   hold off;
   %plot holding current 
   subplot('Position',[0.1,0.10,0.8,0.150]);
   plot(TM, hold_cur, '-r^', 'Markersize', msize, 'Markerfacecolor', 'r') % data in red
   set(gca, 'FontSize', fsize);
 %   v=get(gca, 'YLim');
%    if(v(2) < 0)
%       set(gca, 'Ylim', [v(2) 0]);
%    else
%       set(gca, 'Ylim', [0 v(1)]);
%    end;
    grid;
    ylabel('Holding (nA)');

  
   % plot textual information abstracted from analysis...
   orient landscape
   
   % control printing and closing of window for automatic runs
   % f = 1 creates plot and leaves it up
   % f = 2 creates plot but closes it when done
   % f = 3 creates plot and prints it and then closes it
   if (f == 2)
      print -dljet3;
   end
   if f >= 2  
      close
   end
end

